• Probability is the expression of belief in some future outcome

• A random variable can take on different values with different probabilities

• The sample space of a random variable is the universe of all possible values

• The sample space can be represented by a

  • probability distribution (for discrete variables)
  • probability density function (PDF - for continuous variables)
  • algebra and calculus are used for each respectively
  • the probabilities of an entire sample space always sum to 1.0
• There are many families or forms of distributions or PDFs
  • depends on the nature of the dynamical system they represent
  • the exact instantiation of the form depends on their parameter values
  • we are often interested in statistics in estimating parameters

• Describes the expected outcome of a single event with probability p

• Example of flipping of a fair coin once

\[Pr(X=\text{Head}) = \frac{1}{2} = 0.5 = p \]

\[Pr(X=\text{Tails}) = \frac{1}{2} = 0.5 = 1 - p \]

• If the coin isn't fair then \(p \neq 0.5\)
• However, the probabilities still sum to 1

\[ p + (1-p) = 1 \]

• Same is true for other binary possibilities
  • success or failure
  • yes or no answers
  • choosing an allele from a population based upon allele frequences

• A binomial distribution results from the combination of several independent Bernoulli events

• Example - pretend that you flip 20 fair coins and record the number of heads
• Now repeat that process and record the number of heads for each
• We expect that most of the time we will get approximately 10 heads
• Sometimes we get many fewer heads or many more heads

• The distribution of probabilities for each combination of outcomes is

\[\large f(k) = {n \choose k} p^{k} (1-p)^{n-k}\]

• n is the total number of trials
• k is the number of successes
• p is the probability of success
• q is the probability of not success
• For binomial as with the Bernoulli p = 1-q

• Pretend that you want to create your own binomial distribution
  • you could flip 20 fair coins: 20 independent Bernoulli “trials”
  • you could then replicate this 100 times: 100 “observations”
  • for each one of your replicates you record the number of heads
  • probability of heads (“success”) for any flip of a fair coin is 0.5
• Draw a diagram of what you think the distribution would look like
  • what will the x-axis represent, and what are its values?
  • what will the y-axis represent, and what are its values?
  • where will the 'center of mass' be?
  • what parameter determines the value of the center of mass?
  • what parameter determines the spread of the distribution?

• Using R, simulate these experimental data
  • using the \(rbinom()\) function
  • look at the arguments to see how they map on to the trials, probability of each trial, and number of replicates
  • use the \(hist()\) function to plot the simulated data.
  • what do the y- and x-axes represent?
  • what is the most common outcome, in terms of the number of “successes” per trial? Does this make sense?
• Now just change you script to do the following
  • perform 200 trials for each of the 100 replicates
  • perform 2000 trials for each of the 100 replicate
  • how do the distributions change when you go from 20 to 200 to 2000 trials?

• Note that the binomial function incorporates both the 'and' and 'or' rules of probability
• This part is the probability of each outcome (multiplication)

\[\large p^{k} (1-p)^{n-k}\]

This part (called the binomial coefficient) is the number of different ways each combination of outcomes can be achieved (summation)

\[\large {n \choose k}\] Together they equal the probability of a specified number of successes

\[\large f(k) = {n \choose k} p^{k} (1-p)^{n-k}\]

• Now perform your experiment using R
• Use the rbinom function to replicate what you sketched out for coin flipping
• Be sure to check out the other functions in the binom family
• Make sure you know how you are mapping the parameeters to values
• Take the output and create a histogram
• Does it look similar to what you expected?
• Now change the values around so that
  • You have more or fewer coin flips per trial
  • You replicate the process with more or fewer trials
  • You change the coin from being fair to being slightly biased
  • You change the coin from being slightly to heavily biased

• Another common situation in biology is when each trial is discrete but number of observations of each outcome is observed

• Some examples are
  • counts of snails in several plots of land
  • observations of the firing of a neuron in a unit of time
  • count of genes in a genome binned to units of 500 amino acids in length
• Just like before you have 'successes', but
  • now you count them for each replicate
  • the replicates now are units of area or time
  • the values can now range from 0 to an arbitrarily large number

• For example, you can examine 100 plots of land
  • count the number of snails in each plot
  • what is the probability of observing a plot with 'r' snails is
• Pr(Y=r) is the probability that the number of occurrences of an event y equals a count r in the total number of trials

\[Pr(Y=r) = \frac{e^{-\mu}\mu^r}{r!}\]

• Note that this is a single parameter function because \(\mu = \sigma^2\), and the two together are often just represented by \(\lambda\)

\[Pr(y=r) = \frac{e^{-\lambda}\lambda^r}{r!}\]

• This means that for a variable that is truly Poisson distributed:
  • the mean and variance should be equal to one another, a hypothesis that you can test
  • variables that are approximately Poisson distributed but have a larger variance than mean are often called 'overdispersed'

where \[\large \pi \approx 3.14159\]

\[\large \epsilon \approx 2.71828\]

To write that a variable (v) is distributed as a normal distribution with mean \(\mu\) and variance \(\sigma^2\), we write the following:

\[\large v \sim \mathcal{N} (\mu,\sigma^2)\]

Estimate of the mean from a single sample

\[\Large \bar{x} = \frac{1}{n}\sum_{i=1}^{n}{x_i} \]

Estimate of the variance from a single sample

\[\Large s^2 = \frac{1}{n-1}\sum_{i=1}^{n}{(x_i - \bar{x})^2} \]

The standard deviation is the square root of the variance

\[\Large s = \sqrt{s^2} \]

• Often we want to make variables more comparable to one another
• For example, consider measuring the leg length of mice and of elephants
• Which animal has longer legs in absolute terms?
• Which has longer legs on average proportional to their body size? Which has more variation proportional to their body size?
• A qood way to answer this last question is to use 'z-scores', which are standardized to a mean of 0 and a s.d. of 1
• We can modify any normal distribution to have a mean of 0 and a standard deviation of 1
• Another term for this is the standard normal distribution

\[\huge z_i = \frac{(x_i - \bar{x})}{s}\]

Simulate a population of 10,000 individual values for a variable x:

x <- rnorm(10000, mean=50.5, sd=5.5) 

Take 1000 random samples of size 20, take the mean of each sample, and plot the distribution of these 1000 sample means.

x_sample_means <- NULL
for(i in 1:1000){
x_samp <- sample(x, 20, replace=FALSE)
x_sample_means[i] <- mean(x_samp)
}

For one of your samples, use the equation from the previous slide to transform the values to z-scores.

Plot the distribution of the z-scores, and calculate the mean and the standard deviation.

• Estimation is the process of inferring a population parameter from sample data
• The value of one sample estimate is almost never the same as the population parameter because of random sampling error
• Imagine taking multiple samples, each will be different from the true value
• Most will be close, but some will be far away
• The expected value of a very large number of sample estimates is the value of the parameter being estimated
• Sampling distribution of an estimate
  • all values we might have obtained from our sample
  • probabilities of occurrence
• Standard error of an estimate
  • standard deviation of a sampling distribution
  • measures the precision of the parameter estimate
  • NO ESTIMATE IS USEFUL WITHOUT IT!

x <- c(0.9, 1.2, 1.2, 1.3, 1.4, 1.4, 1.6, 1.6, 2.0, 2.0)

• Use R to make 1000 “pseudo-samples” of size 10 (with replacement), using a for loop as before.
• Name the pseudo-sample object “xboot”, and name the means of the xboot samples “z”.
• Plot the histogram of the resampled means, and calculate the standard deviation of the sample means (the bootstrap SEM) using the sd() function.
• How does it compare with the ordinary standard error of the mean calculated from the original, real sample?

sd(x)/sqrt(10)

• Now take one of the genes from the GacuRNAseq_Subset.csv data and obtain a bootstrapped estimate of the mean expression level.

\[\huge t = \frac{(\bar{y}_1-\bar{y}_2)}{s_{\bar{y}_1-\bar{y}_2}} \]

where

which is the calculation for the standard error of the mean difference

1. Make a dummy data set that contains one continuous and one categorical value (with two levels). Draw individuals of the two different factor levels from normal distributions with slightly different means. Take 100 obsv. per factor level.

2. Now, perform a t-test to see if the two populations are statistically different from one another

boxplot(continuous_variable~cat_variable, dataset name) 
t.test(continuous_variable~cat_variable, dataset name) 

1. Repeat steps 1 and 2 above but use different sample sizes (e.g. 10 , then 100, then 1000, then 10,000 obsv. per factor level)

2. Repeat steps 1 and 2 above but now test between a categorical and continuous variable in the perchlorate data set (or any of your favorite data sets)

• Now, apply the equations for the t-statistic above to your simulated sample to perform the same test, but using a resampling approach
  • hint: go back to your Bootstrapping script
• Use this approach to calculate a distribution of random t statistics assuming the the null hypothesis of no difference is true.
• Now, compare your observed value to your created null distribution.
• What is the probability of that value occurring by chance under the null?
  • This is your p-value! (Assume you’re doing a one-tailed test)
  • hint: Proportion of values in the null distribution equal to or larger than the observed t-value = p-value

• Steps of the randomization test
• First, create a randomized data set in which the measurements are randomly reassigned (without replacement) to the two groups
  • For example, the following
  • SM: 4 0 7 4 2 2 2 1 4 4 0 3 3 4 6 2 4 6 0 0
  • DM: 2 2 3 3 2 3 3 4 1 2 1 4 6 6 2 0
• Second, calculate the test statistic measuring the association between variables
  • difference between group means
  • this is the first bootstrap replicate
• Repeat steps 1 and 2 many times to produce many bootstrap replicates

• Below is the null distribution of the test statistic from 10,000 replications
  • This is produced by the randomized assignment of values to each group (SM or DM)
  • Thus this is the distribution we'd expect under the null hypothesis of no difference
• The observed value from the data is –1.425
• The area in red is the tail beyond this observed value and is therefore the bootstrap p-value

• Mean of x in each case 9 (exact)

• Variance of x in each case 11 (exact)

• Mean of y in each case 7.50 (to 2 decimal places)

• Variance of y in each case 4.122 or 4.127 (to 3 decimal places)

• Correlation between x and y in each case 0.816 (to 3 decimal places)

• Linear regression line in each case \[ y = 3.00 + 0.50x\]

All an be written in the form

response variable = intercept + (explanatory_variables) + random_error

in the general form:

\[ Y=\beta_0 +\beta_1*X_1 + \beta_2*X_2 +... + error\]

where \(\beta_0, \beta_1, \beta_2, ....\) are the parameters of the linear model

Write a script to read in the perchlorate data set. Now, add the code to perform a linear model of two continuous variables. Notice how the output of the linear model is specified to a new variable. Also note that the variables and dataset are placeholders

my_lm <- lm(dataset$variable1 ~ dataset$variable2)

Now look at a summary of the linear model

summary(my_lm)
print(my_lm)

Now let's produce a nice regression plot

plot(dataset$variable1 ~ dataset$variable2, col = “blue”)
abline(my_lm, col = “red”)

Notice that you are adding the fitted line from your linear model Finally, remake this plot in GGPlot2

• zfish_diet_IA.tsv
• Protein - amount of protein in the diet
• Weight - weight of the fish

• SL - standard length of the fish

• Perchlorate_Data.tsv
• Strain - of zebrafish
• Perchlorate Level - treatment of perchlorate
• T4 Hormone Level
• Follicle Area
• Number of Follicles
• Age_Category
• Testes Area

• Testes Stage

Perchlorate_Data <- read.table(‘perchlorate_data.tsv’, header=T, sep=‘/t’)
head(Perchlorate_Data)

x <- Perchlorate_Data$T4_Hormone_Level
y <- Perchlorate_Data$Testes_Area

perc_lm <- lm(y ~ x)
summary (perc_lm)

plot(y ~ x, col = "blue")
abline(perc_lm, col = "red")

Estimation of the variation that is explained by the model (SS_model)

SS_model = SS_total(reduced model) - SS_residual(full model)

The variation that is unexplained by the model (SS_residual)

SS_residual(full model)

• 1 is an outlier for both Y and X
• 2 is not an outlier for either Y or X but has a high residual
• 3 is an outlier in just X - and thus a high residual - and therefore has high influence as measured by Cook's D

• Leverage - a measure of how much of an outlier each point is in x-space (on x-axis) and thus only applies to the predictor variable. (Values > 2*(2/n) for simple regression are cause for concern)

• Residuals - As the residuals are the differences between the observed and predicted values along a vertical plane, they provide a measure of how much of an outlier each point is in y-space (on y-axis). The patterns of residuals against predicted y values (residual plot) are also useful diagnostic tools for investigating linearity and homogeneity of variance assumptions

• Cook’s D statistic is a measure of the influence of each point on the fitted model (estimated slope) and incorporates both leverage and residuals. Values ≥ 1 (or even approaching 1) correspond to highly influential observations.

Let’s look at the zebrafish diet data.

x <- zfish_diet$Protein
y <- zfish_diet$Weight

zfish_lm <- lm(y ~ x)
summary (zfish_lm)

plot(y ~ x, col = "blue")
abline(zfish_lm, col = "red")

hist(residuals(zfish_lm), breaks=30)

plot (residuals(zfish_lm) ~ fitted.values(zfish_lm))
plot (residuals(zfish_lm) ~ x)

Or apply the plot() function to the linear model object directly

plot(zfish_lm)

Figure out what these plots are telling you

How about using the influence.measures function???

influence.measures(zfish_lm)

Do we have any high leverage observations we need to worry about??? Now go back and try to recreate these graphs in GGPlot2

• Residuals - As the residuals are the differences between the observed and predicted values along a vertical plane, they provide a measure of how much of an outlier each point is in y-space (on y-axis). The patterns of residuals against predicted y values (residual plot) are also useful diagnostic tools for investigating linearity and homogeneity of variance assumptions

• Leverage - a measure of how much of an outlier each point is in x-space (on x-axis) and thus only applies to the predictor variable. (Values > 2*(2/n) for simple regression are cause for concern)

• Cook’s D statistic is a measure of the influence of each point on the fitted model (estimated slope) and incorporates both leverage and residuals. Values ≥ 1 (or even approaching 1) correspond to highly influential observations.

• running medians (or less robust running means) generate predicted values that are the medians of the responses in the bands surrounding each observation.
• loess and lowesse (locally weighted scatterplot smoothing) - fit least squares regression lines to successive subsets of the observations weighted according to their distance from the target observation and thus depict changes in the trends throughout the data cloud.
• kernel smoothers - new smoothed y-values are computed as the weighted averages of points within a defined window (bandwidth) or neighbourhood of the original x-values. Hence the bandwidth depends on the scale of the x-axis. Weightings are determined by the type of kernel smoother specified, and for. Nevertheless, the larger the window, the greater the degree of smoothing.
• splines - join together a series of polynomial fits that have been generated after the entire data cloud is split up into a number of smaller windows, the widths of which determine the smoothness of the resulting piecewise polynomial.

x <- zfish$SL
y <- zfish$Weight

size_lm <- lm(y~x)
summary(size_lm)

plot(y~x, col = "blue")
abline(size_lm, col = "red")

Note that when we run the lm() function we use Model I regression. Is this appropriate?

install.packages(“lmodel2”)
library(lmodel2)

size_lm_ModII <- lmodel2(y ~ x)
size_lm_ModII

abline(a=0.1912797, b=0.1097880, lty=2)

You should have generated OLS, MA, and SMA regression models, and the last plots SMA line from parameter estimates

• Say you measured green fluorescent protein abundance in cell culture across several key, carefully controlled temperatures.

• You ultimately want to know whether protein expression changes as a function of temperature in your experiment.

• Read in the gfp_temp.tsv data and perform a regression analysis.

• Address the following questions
  • Which is X and which is Y?
  • Are residual assumptions met?
  • What model of regression should we use in this case and why?

• What is our decision regarding our null hypothesis of no relationship?

• power
• exponential
• polynomial

• Read in the RNAseq_lip.tsv and examine the continuous variables.
• We are interested in whether Gene01 and/or Gene02 expression levels influence lipid content.
• First plot \(Y\) vs \(Gene01\), \(Y\) vs \(Gene02\), then set up and test a multiplicative model.
• What are the parameter estimates of interest in this case?
• What are the outcomes from our hypothesis tests?

RNAseq_Data <- read.table(‘RNAseq_lip.tsv', header=T, sep=‘\t')

y <- RNAseq_Data$Lipid_Conc
g1 <- RNAseq_Data$Gene01
g2 <- RNAseq_Data$Gene02
g3 <- RNAseq_Data$Gene03
g4 <- RNAseq_Data$Gene04

Mult_lm <- lm(y ~ g1*g2)
summary(Mult_lm)

• Now get rid of the interaction term, and set up a purely additive model
• Did any of our estimates change? Why?
• Did the degrees of freedom change? Why?

Add_lm <- lm(y ~ g1+g2)
summary(Add_lm)

• Finally try different combinations of genes

RNAseq_Data <- read.table(‘RNAseq_lip.tsv', header=T, sep=‘\t')

RNAseq_lm_linear <- lm(y ~ g4)
summary (RNAseq_lm_linear)
influence.measures(RNAseq_lm_linear)

RNAseq_lm_poly <- lm(y ~ poly(g4, 2))
summary (RNAseq_lm_poly)

lines(y, fitted(lm(y ~ poly(g4, 2))), type="l", col="blue")
influence.measures(RNAseq_lm_poly)

• Try this for different genes
• Try this for different degree polynomials

• MSresiduals- represents the mean amount of variation unexplained by the model, and therefore the lowest value indicates the best fit.

• Adjusted \(r^2\) - (the proportion of mean amount of variation in response variable explained by the model) is calculated as adj. r2 which is adjusted for both sample size and the number of terms. Larger values indicate better fit.

• Mallow’s Cp - is an index resulting from the comparison of the specific model to a model that contain all the possible terms. Models with the lowest value and/or values closest to their respective p (the number of model terms, including the y-intercept) indicate best fit.

• Akaike Information Criteria (AIC) - there are several different versions of AIC, each of which adds a different constant designed to penalize according to the number of parameters and sample size to a likelihood function to produce a relative measure of the information content of a model. Smaller values indicate more parsimonious models. As a rule of thumb, if the difference between two AIC values (delta AIC) is greater than 2, the lower AIC is a significant improvement in parsimony.

• Schwarz or Bayesian Information Criteria (BIC or SIC) - is outwardly similar to AIC. The constant added to the likelihood function penalizes models with more predictor terms more heavily (and thus select more simple models) than AIC. It is for this reason that BIC is favored by many researchers.

• Again, use the RNAseq_lip.tsv data again.
• Let’s test for an effect of Population on Gene01 expression levels
• First, let’s look at how the data are distributed

RNAseq_Data <- read.table('RNAseq_lip.tsv', header=T, sep='\t')
g1 <- RNAseq_Data$Gene01
Pop <- RNAseq_Data$Population
boxplot(g1~Pop, col=c("blue","green"))

Or, to plot all points:

stripchart(g1~Pop, vertical=T, pch=19, col=c("blue","green"), 
           at=c(1.25,1.75), method="jitter", jitter=0.05)
Pop_Anova <- aov(g1 ~ Pop)
summary(Pop_Anova)

• Below are cases where you are likely to treat factors as random effects
• Whenever your sampling design is nested
  • quadrats within transects
  • transects within woodlots
  • woodlots within districts
• Whenever you divide up plots and apply separate treatments to subplots
• Whenever your replicates are grouped spatially or temporally
  • in blocks
  • in batches
• Whenever you take measurements on related individuals
• Whenever you measure subjects or other sampling units repeatedly

• The F-ratio test for a single-factor ANOVA tests for any difference among groups.
• If we want to understand specific differences, we need further “contrasts”.
• Unplanned comparisons (post hoc):
  • Multiple comparisons carried out after the results are obtained.
  • Used to find where the differences lie (which means differ from which other means)
  • Comparisons require protection for inflated Type 1 error rates:
    • Tukey tests: compare all pairs of means and control for multiple comparisons
    • Scheffé contrasts: compare all combinations of means
• Planned comparisons (a priori):
  • Comparisons between group means that were decided when the experiment was designed (not after the data were in)
  • Must be few in number to avoid inflating Type 1 error rates

• Take the RNAseq data you've examined before and create a new four level genotype by combining genotype and microbiota treatment into a single variable
• Think about how to do this using dplyr functions.

RNAseq_Data <- read.table("RNAseq.tsv", header=T, sep='')

x <- RNAseq_Data$categorical_var
y <- RNAseq_Data$continuous_var1
z <- RNAseq_Data$continuous_var2

• Set up the a priori contrasts specifically testing one group mean against another
• These are just examples - you should figure out the logic of the contrasts

contrasts(x) <- cbind(c(0, 1, 0, -1), c(2, -1, 0, -1), c(-1, -1, 3, -1))

• Confirm that the contrasts are orthogonal

round(crossprod(contrasts(x)), 2)

• Define the contrast labels

rnaseq_data_list <- list(x = list(‘xxx vs. xxx’ = 1, ‘xxx vs. xxx’ = 2, ‘xxx vs. xxx’ = 3))

• Then fit the fixed effect model

RNAseq_aov_fixed <- aov(y ~ x)
plot(RNAseq_aov_fixed)
boxplot(y ~ x)
summary(RNAseq_aov_fixed, split = rnaseq_data_list)

• Remember that this is when you had no hypotheses of differences in means in advance
• Read in the perchlorate data from Week 3
• Let’s assess the effects of the 4 perchlorate levels on T4
• Which perchlorate levels differ in their effect on T4?

perc <- read.table('perchlorate_data.tsv', header=T, sep='\t')

x <- perc$Perchlorate_Level
y <- log10(perc$T4_Hormone_Level)

MyANOVA <- aov(y ~ x)
summary (MyANOVA)
boxplot(y ~ x)

install.packages("multcomp")
library(multcomp)

summary(glht(MyANOVA, linfct = mcp(x = "Tukey")))

• Transmission of the HIV-1 virus via sex workers contributes to the rapid spread of AIDS in Africa
• The spermicide nonoxynol-9 had shown in vitro activity against HIV-1, which motivated a clinical trial by van Damme et al. (2002).
• They tested whether a vaginal gel containing the chemical would reduce the risk of acquiring the disease by female sex workers.
• Data were gathered on a volunteer sample of 765 HIV-free sex-workers in six clinics in Asia and Africa.
• Two gel treatments were assigned randomly to women at each clinic.
• One gel contained nonoxynol-9 and the other contained a placebo.
• Neither the subjects nor the researchers making observations at the clinics knew who received the treatment and who got the placebo.

The goal of experimental design is to eliminate bias and to reduce sampling error when estimating and testing effects of one variable on another.

• To reduce bias, the experiment included:
  • Simultaneous control group: study included both the treatment of interest and a control group (the women receiving the placebo).
  • Randomization: treatments were randomly assigned to women at each clinic.
  • Blinding: neither the subjects nor the clinicians knew which women were assigned which treatment.
• To reduce the effects of sampling error, the experiment included:
  • Replication: study was carried out on multiple independent subjects.
  • Balance: number of women was nearly equal in the two groups at every clinic.
  • Blocking: subjects were grouped according to the clinic they attended, yielding multiple repetitions of the same experiment in different settings (“blocks”).

• In clinical trials either a placebo or the currently accepted treatment should be provided.
• In experiments requiring intrusive methods to administer treatment, such as
  • injections
  • surgery
  • restraint
  • confinement
• the control subjects should be perturbed in the same way as the other subjects, except for the treatment itself, as far as ethical considerations permit.
• The “sham operation”, in which surgery is carried out without the experimental treatment itself, is an example.
• In field experiments, applying a treatment of interest may physically disturb the plots receiving it and the surrounding areas, perhaps by trampling the ground by the researchers.
• Ideally, the same disturbance should be applied to the control plots.

• In a single-blind experiment, the subjects are unaware of the treatment that they have been assigned.
• Treatments must be indistinguishable to subjects, which prevents them from responding differently according to knowledge of treatment.
• Blinding can also be a concern in non-human studies where animals respond to stimuli
• In a double-blind experiment the researchers administering the treatments and measuring the response are also unaware of which subjects are receiving which treatments
  • Researchers sometimes have pet hypotheses, and they might treat experimental subjects in different ways depending on their hopes for the outcome
  • Many response variables are difficult to measure and require some subjective interpretation, which makes the results prone to a bias
  • Researchers are naturally more interested in the treated subjects than the control subjects, and this increased attention can itself result in improved response

• Replication is the assignment of each treatment to multiple, independent experimental units.
• Without replication, we would not know whether response differences were due to the treatments or just chance differences between the treatments caused by other factors.
• Studies that use more units (i.e. that have larger sample sizes) will have smaller standard errors and a higher probability of getting the correct answer from a hypothesis test.
• Larger samples mean more information, and more information means better estimates and more powerful tests.
• Replication is not about the number of plants or animals used, but the number of independent units in the experiment. An “experimental unit” is the independent unit to which treatments are assigned.
• The figure shows three experimental designs used to compare plant growth under two temperature treatments (indicated by the shading of the pots). The first two designs are un-replicated.

• A study design is balanced if all treatments have the same sample size.
• Conversely, a design is unbalanced if there are unequal sample sizes between treatments.
• Balance is a second way to reduce the influence of sampling error on estimation and hypothesis testing.
• To appreciate this, look again at the equation for the standard error of the difference between two treatment means.

• For a fixed total number of experimental units, n1 + n2, the standard error is smallest when n1 and n2 are equal.
• Balance has other benefits. For example, ANOVA is more robust to departures from the assumption of equal variances when designs are balanced or nearly so.

• Nested ANOVA or nested design
  • factors might be hierarchical - in other words nested - within one another
  • The sources of variance are therefore hierarchical too
• The factorial ANOVA design is the most common experimental design used to investigate more than one treatment variable
  • In a factorial design every combination of treatments from two (or more) treatment variables is investigated.
  • The main purpose of a factorial design is to evaluate possible interactions between variables.
  • An interaction between two explanatory variables means that the effect of one variable on the response depends on the state of a second variable.

• Nested designs are hierarchical
  • often contain sub-replicates that are random, uncontrolled, nuisance effects
  • but the nested factors can be of interest too
• Factorial designs are
  • all pairwise combinations,
  • and often involve all combinations of factor levels
  • when each factor is fixed interactions can be assessed
• Completely nested designs therefore have no interaction terms, whereas factorial designs do
• Mixed models can have a combination of fixed and random factors that are more complicated

• Example 1: Study of “repeatability” (simple nested design)
• The walking stick, Timema cristinae, is a wingless herbivorous insect on plants in chaparral habitats of California.
• Nosil and Crespi (2006) measured individuals using digital photographs.
• To evaluate measurement repeatability they took two separate photographs of each specimen.
• After measuring traits on one set of photographs, they repeated the measurements on the second set.

andrew_data <- read.table('andrew.tsv', header=T, sep=‘\t')
head(andrew_data)

• There are four variables: ‘TREAT’, ‘PATCH’, ‘QUAD’ and ‘ALGAE’
• The main effect factor is TREAT
• Make a simplified factor called TREAT2, in which 0% and 33% are a level called “low” and 66% and 100% are “high”

andrew_data$TREAT2 <- factor(c(rep(“low”,40),rep(“high”,40))

• The nested factor is PATCH - also need to turn this into a factor

andrew_data$PATCH <- factor(andrew_data$PATCH)

• In this case, our response variable is ALGAE
• Look at the distribution of ALGAE for the two levels of TREAT2 using boxplots based on the patch means, which are the replicates in this case.

andrew.agg <- with(andrew_data, aggregate(data.frame(ALGAE), 
                  by = list(TREAT2=TREAT2, PATCH=PATCH), mean)

library(nlme)
andrew.agg <- gsummary(andrew_data, groups=andrew_data$PATCH)

boxplot(ALGAE ~ TREAT2, andrew.agg)

• Evaluate assumptions based on the boxplots
• Is the design balanced (equal numbers of sub-replicates per PATCH)?

• Run the nested ANOVA:

nested.aov <- aov(ALGAE ~ TREAT2 + Error(PATCH), data=andrew_data)
summary(nested.aov)

• Do we detect an effect of TREAT2 (high vs low sea urchin density)?
  
• Estimate variance components to assess relative contributions of the random factors

library(nlme)
VarCorr(lme(ALGAE ~ 1, random = ~1 | TREAT2/PATCH, andrew_data))

• Calculate the % of variation due to between-treatment differences vs. due to among patches within treatment differences.
• See pg. 302 in Logan if you need help.
• What do these variance component estimates tell us???

• For example, Relyae (2003) looked at how a moderate dose (1.6mg/L) of a commonly used pesticide, carbaryl (Sevin), affected bullfrog tadpole survival.
• In particular, the experiment asked how the effect of carbaryl depended on whether a native predator, the red-spotted newt, was also present.
• The newt was caged and could cause no direct harm, but it emitted visual and chemical cues to other tadpoles
• The experiment was carried out in 10-L tubs (experimental units), each containing 10 tadpoles.
• The four combinations of pesticide treatment (carbaryl vs. water only) and predator treatment (present or absent) were randomly assigned to tubs.
• The results showed that survival was high except when pesticide was applied together with the predator.
• Thus, the two treatments, predation and pesticide, seem to have interacted.

continuous response and two main effect variables

rnadata <- read.table('RNAseq.tsv', header=T, sep='')
gene <- rnadata$Gene80
microbiota <- rnadata$Microbiota
genotype <- rnadata$Genotype

make new “pseudo factor,” combining genotype and microbiota

gxm <- interaction(genotype,microbiota)
levels(gxm)
boxplot(gene ~ gxm)

specify the following 2 contrasts

contrasts(gxm) <- cbind(c(2, -1, 0, -1), c(-1, -1, 3, -1))

Fit the factorial linear model

rna_aov <- aov(gene ~ gxm)

Examine the ANOVA table, using supplied contrasts. Figure out the appropriate titles to give them.

summary(rna_aov, split = list(gxm = list('xxx'=1,'xxx'=2)))

What does the contrast summary tell you about the nature of the interaction?

rnadata <- read.table('RNAseq.tsv', header=T, sep='')
head(rnadata)

variables excluding first 5 and last 5 observations

gene <- rnadata$Gene80[6:75] 
microbiota <- rnadata$Microbiota[6:75]
genotype <- rnadata$Genotype[6:75]
boxplot(gene ~ microbiota)
boxplot(gene ~ genotype)
boxplot(gene ~ microbiota*genotype)

Estimate the variance components using Restricted Maximum Likelihood (REML)

library(lme4)
lmer(gene ~ 1 + (1 | microbiota) + (1 | genotype) + (1 | microbiota:genotype))

Based on the REML sd estimates, what are the relative contributions of the factors to total variance in gene expression?

longevity_data <- read.table('longevity.csv', header=T, sep=',')
head(longevity_data)

Variables

long <- longevity_data$LONGEVITY
treat <- longevity_data$TREATMENT
thorax <- longevity_data$THORAX

• check to see if the covariate should be included

boxplot(long ~ treat)
plot(long ~ thorax)

• assess assumptions of normality and homogeneity of variance

plot(aov(long ~ thorax + treat ), which = 1)

• †ry it again with a transformed response variable

plot(aov(log10(long) ~ thorax + treat ), which = 1)

• visually assess linearity, homogenetiy of slopes and covariate range equality

library(lattice)
print(xyplot(log10(long) ~ thorax | treat, type = c("r", "p")))

• formally test homogenetiy of slopes by testing the interaction term

anova(aov(log10(long) ~ thorax*treat))

• formally test covariate range disparity by modeling the effect of the treatments on the covariate

anova(aov(thorax ~ treat))

• FINALLY, set up contrasts, fit the additive model and visualize the results (pg. 459 and 460 of your Logan book)
• Summarize the trends in a nice plot (pg. 461 of your Logan book)

1. Distributions of data
   • location
   • spread
   • shape
2. Associations between variables
   • relationship among two or more variables
   • differences among groups in their distributions

1. Distributions of data
   • bar graph
   • histogram
   • box plot
2. Associations between variables
   • pie chart
   • grouped bar graph
   • mosaic plot
   • box plot
   • scatter plot
   • dot plot 'stripchart'